Developer CD Series 1994…tember: Reference Library/Dev.CD Sep 94.toast/Periodicals/develop/develop Issue 11/develop 11 code/The NetWork Project/Examples (Sources)/NetSim/HistogramUnit.inc

Text File | 1992-07-15 | 7.9 KB | 315 lines | [TEXT/MPS ]

    USES sane;{ HistogramUnit.inc     © Copyright G. Sawitzki, 1988-1991}
    {$S Histogram}

PROCEDURE initStat(VAR stat:tStatType;name:str15);
BEGIN
    WITH stat DO BEGIN
        count:=0;
        mean:=0;
        ssq:=0;
        min:=+Inf; max:=-inf;
        id:=name;
    END;
END;

{Subroutine  AddStat 1.2}
{  For evaluation of mean and sum of deviation squares}
{ this procedure executes one single iteration}
{ obeying a provisional mean algorithm.}
{ At "w" the actual value is expected,}
{ at "stat" mean,SSQ and number of the preceding ones}
{ according the statlab-convention.}
{ type stattyp = record}
{        count : integer;}
{        mean  : extended;}
{        ssq   : extended}
{                end;}
{}
{Updating history:}
{1.2 added mi/max. gs.}


PROCEDURE addstat;
    {(w : extended;    var    stat : stattyp    );    }
VAR
    xDiff: extended;
BEGIN
    IF (ClassExtended(w) in [SNAN,QNAN,Infinite]) THEN  SysBreakStr('w');
    WITH stat DO
    BEGIN
        count := count + 1;
        xDiff := w - mean;
        mean := mean + xDiff / count;
    IF (ClassExtended(mean) in [SNAN,QNAN,Infinite]) THEN  SysBreakStr('mean');
        ssq := ssq + xDiff * (w - mean);
        IF w > max THEN
        BEGIN
            max := w;
        END;
        IF w < min THEN
        BEGIN
            min := w;
        END;
    END;
END;

PROCEDURE addbistat;

{ (x, y : extended;var xstat, ystat : tStatType;    var    xyssq : extended    ) ;   }


    {Fuction: accumulate (x,y) into the statistics}
    {Input: (x,y): new data}
    {       xstat,ystat,xyssq: current statistics}
    {Output:  xstat,ystat,xyssq: updated statistics}


    {Updating by method of provisional means}
    {BMDP Statistical Software, Appendix A.2}
    {xyssq: sum of (x-xmean)(y-ymean)}
    {no checks performed, count overflow will not be detected}

VAR
    xDiff, ydiff: extended;
BEGIN
    {update statistics for x}
    WITH xstat DO
    BEGIN
        count := succ(count);
        xDiff := x - mean;
        mean := mean + xDiff / count;
        ssq := ssq + xDiff * (x - mean);
    END;
    {update statistics for y and cross-products}
    WITH ystat DO
    BEGIN
        count := count + 1;
        yDiff := y - mean;
        mean := mean + yDiff / count;
        ssq := ssq + yDiff * (y - mean);
        xyssq := xyssq + xdiff * (y - mean);
    END;
END;


PROCEDURE initHistogram (VAR hist: histtype;name:str15);
    {(   var    hist : histtype   );}
VAR
    I: integer;
BEGIN

    WITH hist DO
    BEGIN
        min := +INF;
        max := -Inf;
        BinWidth := Inf;
        count := 0;
        maxbincount := 1;{should be 1 to allow scaling for empty histogram}
        FOR I := 0 TO maxclass DO
        cnt[I] := 0;
        id:=name;
    END;
END;

FUNCTION histogramcell(val: extended; VAR hist: histtype):integer;
VAR tempcell:integer;templim,delta:extended;
BEGIN
    WITH hist DO
    BEGIN
        IF val>max THEN histogramcell:=maxint
        ELSE IF val<=min THEN histogramcell:=-1
        ELSE BEGIN
         {$IFC PARANOIA} IF binwidth<=0 THEN debugstr('bin widht <=0'); {$ENDC}
        {$IFC PARANOIA} IF binwidth=Inf THEN debugstr('bin widht Inf'); {$ENDC}
            tempcell:=trunc((val-min)/binwidth);
            templim:=min+tempcell*binwidth;        {left boundary of this cell}
            delta:=val-min+tempcell*binwidth;    {?? sometimes val<templim}
            IF val<=templim THEN
            tempcell:=tempcell-1 {we checked for min, so tempcell>0}
            ELSE
            BEGIN
                templim:=templim+binwidth;
                IF val>templim THEN Debugstr('larger')
            END;
      {$IFC PARANOIA} IF tempcell<0 THEN debugstr('tempcell negative'); {$ENDC}

            histogramcell:=tempcell;
        END;
    END;{with}
END;

PROCEDURE addhist (val: extended; VAR hist: histtype);

PROCEDURE secondvalue (val: extended);
VAR
    countmax, countmin: longint;
    slop,valmin, valmax: extended;
    cell: integer;
BEGIN {second value}
    WITH hist DO
    BEGIN
        IF val > max THEN
        BEGIN
            min := max;
            countmin := count - 1;    {we already incremented count for val}
            max := val;
            countmax := 1;
        END
        ELSE
        BEGIN
            countmax := count - 1;
            min := val;
            countmin := 1;
        END;

        {min and max are defined & correct}
        valmin := min;                {temporary - copy from misused histogram entries}
        valmax := max;

        {find initial scale. try with quarter range slop}
        slop:=(max - min) / 4;
        min := valmin - slop;
        IF (min < 0) & (valmin >0) THEN min := 0;    {fix on zero}
        max := valmax + slop ;
        IF (max > 0) & (valmax < 0) THEN max := 0;    {fix on zero}

        {$IFC PARANOIA} IF min>max THEN debugstr('bad max/min'); {$ENDC}

        binwidth := (max - min) / (maxClass + 1);

        {scales are set up. Now store everything to correct cell}
        cell := histogramcell(valmin,hist);
        cnt[cell] := countmin;
        IF countmin > maxbincount THEN
        maxbincount := countmin;

        cell := histogramcell(valmax,hist);
        cnt[cell] := countmax;
        IF countmax > maxbincount THEN
        maxbincount := countmax;

    END;
END;{second value}

PROCEDURE movebins(frombin,tobin,nrbins:integer);
    {shift bins, and adjust limits. should work in both directions.
    All bins outside the moved range are considered invalid and
    the bin counts will be set to zero. On the Macintosh, we could
    use CMove to speed up. For portability, we use a loop.}
VAR bin,delta:integer;
BEGIN
    IF nrbins<=0 THEN
    {error. should not occur. NOOP for now}
    ELSE WITH hist DO BEGIN
        delta:=tobin-frombin;
        IF frombin<tobin THEN {shift right. delta >0}
        BEGIN
            FOR bin:=frombin+nrbins-1 DOWNTO frombin DO
            cnt[bin+delta]:=cnt[bin];
        END ELSE IF frombin>tobin THEN {shift left. delta <0}
        BEGIN
            FOR bin:=frombin TO frombin+nrbins-1 DO
            cnt[bin+delta]:=cnt[bin];
        END;
        max:=max-delta*binwidth;
        min:=min-delta*binwidth;
        {set remaining bin counts to zero}
        FOR bin:=0 TO tobin-1 DO cnt[bin]:=0;
        FOR bin:=tobin+nrbins TO maxclass DO cnt[bin]:=0;
    END;
END;

PROCEDURE newcell (val: extended);
VAR
    cell, fromclass,toclass: integer;
    usedcell,delta: integer;
    {val not in range min<val≤max}
BEGIN
    WITH hist DO
    BEGIN
        IF val<= min THEN
        BEGIN
            {left needed. try a shift first}
            usedcell:=maxclass;
            WHILE cnt[usedcell]=0 DO usedcell:=usedcell-1;
            IF usedcell<maxclass THEN BEGIN    {a free cell at right}
                delta:=(maxclass-usedcell) DIV 2;
                IF delta=0 THEN delta:=1;
                MoveBins(0,delta,usedcell+1);
            END
            ELSE BEGIN{left needed. change scale}
                fromclass := maxclass;
                FOR cell := maxclass DOWNTO (maxclass DIV 2) + 1 DO
                BEGIN
                    cnt[cell] := cnt[fromclass] + cnt[fromclass - 1];
                    fromclass := fromclass - 2;
                    IF cnt[cell] > maxbincount THEN
                    maxbincount := cnt[cell];
                END;{shift right}
                FOR cell := (maxclass DIV 2) DOWNTO 0 DO
                cnt[cell] := 0;
                binwidth := binwidth * 2;
                min := min - ((maxclass DIV 2) + 1) * binwidth;
            END {left needed. change scale}
        END
        ELSE
        BEGIN {right needed. try a shift first}
            usedcell:=0;
            WHILE cnt[usedcell]=0 DO usedcell:=usedcell+1;
            IF usedcell> 0 THEN BEGIN    {a free cell at left}
                delta:=usedcell DIV 2;
                IF delta=0 THEN delta:=1;
                movebins(usedcell,usedcell-delta,maxclass-usedcell+1);
            END ELSE
            BEGIN  {right needed. change scale}
                fromclass := 0;
                FOR cell := 0 TO (maxclass DIV 2) DO
                BEGIN
                    cnt[cell] := cnt[fromclass] + cnt[fromclass + 1];
                    fromclass := fromclass + 2;
                    IF cnt[cell] > maxbincount THEN
                    maxbincount := cnt[cell];
                END;{shift left}
                FOR cell := (maxclass DIV 2) + 1 TO maxclass DO
                cnt[cell] := 0;
                binwidth := binwidth * 2;
                {change scale}
                max := min + binwidth * (maxclass + 1);
            END; {right needed. change scale}
        END; {right needed.}
    END;{with hist}
END;

(* main part procedure addHist *)
VAR
    cell: integer;
    countmin, countmax: longint;
BEGIN
    WITH hist DO
    BEGIN
        count := count + 1;

        IF binwidth < Inf THEN    {ordinary case: binwidth defined}
        BEGIN
            WHILE (val <= min) OR (val > max) DO newcell(val);
            cell := histogramcell(val,hist);
{$IFC paranoia} IF (cell<0) OR (cell>maxclass) THEN debugstr('bad histogramcell'); {$ENDC}
            cnt[cell] := cnt[cell] + 1;
            IF cnt[cell] > maxbincount THEN
            maxbincount := cnt[cell];
        END    {ordinary case: binwidth defined}
        ELSE
        BEGIN {exception case: binwidth undefined}
       IF val = max THEN {tie on first value: noop (we already increased count}
            ELSE
            BEGIN
                IF max = -inf THEN {first observation}
                max := val
                ELSE
                secondValue(val);
            END;
        END;
    END;{with hist}
END; (* addhist *)